Time Series Forecasting

Problem Statement

Neo Bank of India (a fictional bank) receives a lot of requests for its various finance offerings that include housing loans, two-wheeler loans, real estate financing, and microloans. The number of applications received is something that varies a lot with the season. Going through these applications is a manual and tedious process.

Accurately forecasting the number of cases received can help with resource and manpower management resulting in quick response on applications and more efficient processing. You have been appointed with the task of forecasting daily cases for the next 1 month for 2 different business segments aggregated at the country level keeping in consideration the following major Indian festivals (inclusive but not exhaustive list): Diwali, Dussehra, Ganesh Chaturthi, Navratri, Holi, etc.

Data

The dataset is divided into 2 components, namely train data and test data. The training data consists of dependent variables such as date of application (application_date), business segment (segment 1 or 2), anonymized id for a branch at which application was received (branch_id), State in which application was received (state), zone of state in which application was received (zone). The target variable here is the number of cases/applications received (no_of_applicants).

The test data consists of a unique id for each sample in the test set (id), Date of application (application_date), business segment (segment).

Sample submission consists of id with the number of applicants forecasted for next month for business segments 1 and 2.

Evaluation Metrics

Evaluation for the given case is MAPE - Mean Absolute Percentage Error

Data Preprocessing

In the given data as a part of pre-processing we divide the data based on the business segments to analyze them individually for the missing values we use Interpolation to fill in the missing values.We combined the dataset with Holiday dataset downloaded from Kaggle through left join. This final dataset is used for further analysis. The Holiday dataset consists of all the Indian holidays for a period of 2017 to 2019.The data type of the Application date is also modified to ‘datetime’ type and the dataframes for indovidual segments are converted to ‘tsibble’ object for time series analysis.

Exploratory Data Analysis

Overall we observe that the data from the period 2017 to 2019 tends to be seasonal with peaks during the mid-months.The number of applicants in 2017 and early 2019 tend to be comparatively low and inconsistent with the trend observed during the late-2017 and the year 2018, this may be due to the recovery from the demonetization in November-2016.This inconsistent data may lead to noise while forecasting. On dividing based on segments

We observe that majority of the applications belong to Segment 1 and the zone to receive highest number of applications is the East zone. It is observed that the state of Maharashtra has highest number of application across is branches from 2017-2019.

segment_count_plot 

zone_count_plot

data_grpp<-train %>% group_by(application_date,segment) %>% 
       summarise(no_of_applicants_sum  = mean(no_of_applicants))

`summarise()` has grouped output by 'application_date'. You can override using the `.groups` argument.

View(data_grpp)

#Filter data frames based on segment, then convert them into tsibble for segment 1
data_11<-data_grpp %>% filter(segment %in% c(1))
data_11<-data_11 %>%mutate(application_date = as.Date(application_date)) %>%
  as_tsibble(index = application_date)

Checking for seasonality, trend and cycles

The below graph shows the number of application against the time period between 2017-19 under segment 1.We notice that the data tends to be stationary with outliers in the region between 2017-2018.Though this a visual conclusion the hypothesis can be validated by testing for stationarity using the ADF(Augmented Dicker-Fuller)Test. Between the period 2017-2018 we observe a linear trend in number of applications received.The trend doesn’t persist for a long time since there is a sudden deceleration in the number of applications during April to May-2018. The applications received follows a cyclical pattern which are usually generated by the autoregressive structure of the series.

Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

data_1%>%
  gg_season(no_of_applicants_sum, period="year") +
  labs(y = "Total number of applications",
       title = "Yearly Seasonal plot: Number of applications in Neo bank")

Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

Compared to the Segment 1 , the number of applications received under Segment 2 are most consistent over the given periods and it follows strict cyclical pattern.We observe a monthly seasonality with highest number of applicants during the middle of the months.

 autoplot(data_2, no_of_applicants_sum) +
   labs(y = "Total number of applications",
       title = "Number of applications in Neo bank for Segment 2")

Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

# b) Yearly seasonality

train_ts %>%
  gg_season(no_of_applicants, period="year") +
  labs(y = "Total number of applications",
       title = "Yearly Seasonal plot: Number of applications in Neo bank")

data_2%>%
  gg_season(no_of_applicants_sum, period="year") +
  labs(y = "Total number of applications",
       title = "Yearly Seasonal plot: Number of applications in Neo bank")

Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

The above graph shows the cyclical nature of the data across the months for 2017 to 2019.

# # 3. Seasonal sub series plots
# 
# options(repr.plot.width=57 ,repr.plot.height= 10)
# data_1 %>%
#   gg_subseries(no_of_applicants_sum) +
#   labs(
#     y = "Total number of applications",
#     title = "Number of applications in Neo bank"
#   )

#Segment+ Zone 
ggplot(Segment1, aes(x = application_date, y = no_of_applicants)) +
geom_line() +
facet_wrap( ~ segment + zone, scales = "free_y", ncol = 5)+
theme(axis.text.x=element_text(angle=90))

From the data we observe the key region for maximum no.of applications North for the future where the bank may expect higher applications.

# data_2 %>%
#   gg_subseries(no_of_applicants_sum) +
#   labs(
#     y = "Total number of applications",
#     title = "Number of applications in Neo bank"
#   )

month_summary <-Segment1 %>% 
                    group_by(Month) %>%
                    summarise(no_of_applicants= sum(no_of_applicants),.groups='drop') 
month_summary %>%
         ggplot(aes(x = reorder(Month, no_of_applicants), y = no_of_applicants, fill = Month)) +
         geom_col() +
         scale_fill_hp_d(option = "RonWeasley") +
         scale_fill_brewer(palette = "Set3") +
         scale_y_continuous(limits = c(0, 300000), expand = c(0,0)) +
         labs(title = "Total number of applicants in  each month", x = "Month", y = "Total applicants") +
         theme_classic() 

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

We observe that the bank receives the highest number of applications in October across the three years. Though may vary region wise and state wise.We observe maximum banking activity starting from end of the financial year.

data%>%
  pivot_wider(values_from=no_of_applicants, names_from=state)%>% 
  GGally::ggpairs(columns = 2:9)

5. Lag plots

COMPONENTS OF THE TIME SERIES FOR EACH SEGMENT

We attempt to analyze three main components of the time series to analyze how they behave over time-

Seasonal Component

Trend Cycle Component

Remainder Component

SEGMENT1

For Segment 1 we observe a linear trend from the time period 2017 to 2019.There is an inconsistent yearly seasonality since the gaps are wide at certain points.Throughout the year we see outliers mainly in the Quarter 3.The yearly seasonality shows that the peak business period to be mainly in the months of August-September-October-November for Segment 1 customers.The residual part removes the residual and trend component from the data and seasonally adjusted data shows a variation in the number of application received thus the seasonality plays a major role in the business problem under consideration for Segment 1 customers.

SEGMENT 2

components(dcmp) %>% autoplot()

On decomposing the data under Segment 2 we observe that there is a monthly seasonality we observe through out the years. The weekly seasonality is also significant since the gaps are neither close by nor wide. The trend that we observe is mostly horizontal through each year and tends to have a long term upward trend from 2017 through 2018 from 2018 we observe a downward trend.On removing the seasonality and trend the seasonally adjusted data shows more number of changing levels/points.

CHOICE OF MODELS

ACF and PACF plots

If the autocorrelation plot of acf crosses the dashed blue line, it means that specific lag is significantly correlated with current series.

From here, we can observe that there is approx peak at 0,30,60,90… lags. The maximum at lag 0,30,60.. indicates a positive relationship with monthly cycle.

Here, instead of tailing-off if there is a slow decay we can say it to be a non- stationary series. Examining ACF and PACF helps us to estimate plausible models and select p,d and q. Differencing may be needed if there is slow- decay in the ACF. Several models may need to be estimated. So, here we will try out different combinations of p,d and q to get model with lower AIC or BIC score. AIC and BIC are two measures of goodness of fit. They measure the trade-off between model fit and complexity of data. A lower AIC or BIC value indicates a better fit.

#install.packages('patchwork')
library(patchwork)

a+ b +c +d

#TSLM Modeling

TSLM modeling: In this case the model function is TSLM() (time series linear model), the response variable is no_of_applicants and it is being modeled using trend() (a “special” function specifying a linear trend when it is used within TSLM()). First of all linear model on segment 1 is applied.

Here, we can observe approx value of future predictions for next 30 days.

Next, same process is applied on segment 2 dataset, and approx number of applicants is observed.

From here, we can observe approx values of number of applicants for next 30 days for segment 2. Each row corresponds to one forecast period for each day. The no_of_applicants column contains the forecast distribution, while the .mean column contains the point forecast. The point forecast is the mean (or average) of the forecast distribution.

Warning in file(file, ifelse(append, "a", "w")) :
  cannot open file 'C:\Users\asus\Desktop\College\3rd semester\Time series-Akhileshwar\tslm_final1.csv': No such file or directory
Error in file(file, ifelse(append, "a", "w")) : 
  cannot open the connection

Finally, we combine dataset for next 30 days for segment 1 and 2.

Dickey fuller test

Error in adf.test(data_1$no_of_applicants_sum) : 
  could not find function "adf.test"

Here, null hypothesis is that time series is not stationary.If null hypothesis is rejected, then target variable is stationary. If it’s not rejected, we should do differencing and then perform dickey-fuller test again. Here, p-value is less than 0.05, so we can reject null hypothesis and say that series is stationary enough to do any kind of time series modeling.

#AR In order to fit the AR model, we have initially converted the number of applicants column to a vector. To see the lag value to be used in the AR model, we have used the above function that has given a value of 7. Hence, we fit a AR(7) model.

Using the predict function, the total number of applications for the next 30 days is forecasted. The (7,0,0) component in arima function is same as fitting a 7th order Autoregressive model. The AIC score for segment 1 data is equal to

AIC(seg1_ar)

[1] 7307.939

The steps are similarly followed for segment 2 data as well. The AIC score for Segment 2 is found to be

[1] 11576.89

Arima Modeling

While exponential smoothing models are based on a description of the trend and seasonality in the data, ARIMA models aim to describe the autocorrelations in the data. White noise describes the assumption that each element in a series is a random draw from a population with zero mean and constant variance. Autoregressive(AR) and moving average(MA) models correct for violations of this white noise assumption.

As our model satisfies stationarity, we move on to applying arima modeling on it. ## Applying arima modeling on segment 1

Series: no_of_applicants_sum 
Model: ARIMA(0,1,3)(0,0,2)[7] 

Coefficients:
          ma1      ma2      ma3    sma1    sma2
      -0.7880  -0.1167  -0.0560  0.0630  0.0853
s.e.   0.0353   0.0466   0.0362  0.0369  0.0332

sigma^2 estimated as 548.5:  log likelihood=-3633.84
AIC=7279.69   AICc=7279.79   BIC=7307.76

Applying arima modeling on segment 2

Series: no_of_applicants_sum 
Model: ARIMA(0,1,0)(0,0,1)[7] 

Coefficients:
        sma1
      0.3526
s.e.  0.0341

sigma^2 estimated as 83552:  log likelihood=-5760.52
AIC=11525.04   AICc=11525.05   BIC=11534.44

Here, we observe that we get better predictions from ARIMA than we got from TSLM model. This can be observed from plots as well.

ARIMA (Without seasonality)

Though we got our data as stationary, but plots appeared to have some sort of seasonality. So, we went ahead to use arima by removing seasonality. In segment 1

Series: no_of_applicants_sum 
Model: ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.1892  -0.9535
s.e.  0.0384   0.0142

sigma^2 estimated as 513.8:  log likelihood=-3609.18
AIC=7224.35   AICc=7224.38   BIC=7238.39

Here, we observe that AIC and BIC value decreases, which shows that this is a better model if we remove seasonality.

Removing seasonality from segment 2

Series: no_of_applicants_sum 
Model: ARIMA(1,1,1)(0,0,2)[7] 

Coefficients:
          ar1     ma1    sma1     sma2
      -0.0446  0.0823  0.0615  -0.2470
s.e.   0.4993  0.4949  0.0330   0.0301

sigma^2 estimated as 63097:  log likelihood=-5644.87
AIC=11299.74   AICc=11299.82   BIC=11323.24

Here, we observe that AIC and BIC value both decreases when we remove seasonality from data, which shows that this is a better model.

Checking for Trend in data

First off, I’d take the difference of time series to de-trend it, it has a considerable stochastic trend. The most obvious sign of this, beside the time series having a steady upward rise, is that the ACF components takes a long time to die out. While we can fit a model to the data as is with a drift term, it’s much easier to read the ACF and PACF plots when the trend isn’t there. So, we also tried to overcome presence of small amount of trend by two methods(differencing and removing trend component from stl decomposition as well).

Data after differencing

Here, we try to observe ACF and PACF model after differencing.

Warning: Removed 1 row(s) containing missing values (geom_path).
Warning: Removed 1 rows containing missing values (geom_point).

Further, we just work on remainder time series, by removing trend and seasonality.

Removing trend and seasonality from data

We start with segment 1

Series: no_of_applicants_sum 
Model: ARIMA(2,0,0) 

Coefficients:
          ar1      ar2
      -0.1770  -0.2877
s.e.   0.0339   0.0339

sigma^2 estimated as 318.6:  log likelihood=-3422.64
AIC=6851.28   AICc=6851.31   BIC=6865.32

For segment 2

Series: no_of_applicants_sum 
Model: ARIMA(0,0,2)(0,0,2)[7] 

Coefficients:
         ma1      ma2    sma1     sma2
      0.2541  -0.1037  0.0741  -0.2111
s.e.  0.0472   0.0539  0.0340   0.0320

sigma^2 estimated as 29586:  log likelihood=-5343.5
AIC=10697.01   AICc=10697.08   BIC=10720.52

Removing trend and seasonality, aic and bic value reduces a lot, but when we observe it’s predictions we see that it has negative values, which doesn’t make sense. So, instead we should try INARIMA model here, to make sure that our output is positive only which can be done as future work.

HOLT WINTER’S EXPONENTIAL SMOOTHING

This method uses two parameters, one for the overall smoothing and the other for the trend smoothing equation. The method is also called double exponential smoothing or trend-enhanced exponential smoothing.For the given data we observe the data the component are additive hence we follow the additive method. This method is used since the forecast is for a short period and it doesn’t work well for long term predictions.

Here we observe the error rate increasing over time and after differentiation too the model captures only

library(patchwork)
e=autoplot(data_1.dif)
f=autoplot(ses.data_1.dif)
e

f

SARIMAX

Seasonal Auto-Regressive Integrated Moving Average with exogenous factors, or SARIMAX, is an extension of the ARIMA class of models. SARIMAX has the capability to handle exogenous variable like holidays. We can also get your own domain-specific features if you need to. In our case, we combined the original data with Indian holidays and built a SARIMAX model.

Initially the binary variables in Holiday column was converted to 1s and 0s from yes and no. This data is added as the exogenous variable into the ARIMAX model, making it SARIMAX.We use the method Maximum Likelihood Estimate to make the predictions. Fitting SARIMAX model with seasonal ARIMA order (2,2,0) for segment 1 and (4,2,0) for segment 2.

For Segment 1:

model_ts

Call:
arimax(x = (inf), order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12), 
    xreg = rinf, method = "ML")

Coefficients:
          ma1     sma1      xreg
      -0.9268  -1.0000  -14.5327
s.e.   0.0172   0.0163    6.5784

sigma^2 estimated as 604.9:  log likelihood = -3261.34,  aic = 6528.68

And a BIC score of :

BIC(model_ts)

[1] 6548.88

For Segment 2

model_ts1

Call:
arimax(x = (iinf), order = c(4, 2, 0), seasonal = list(order = c(4, 2, 0), period = 12), 
    xreg = rrinf, method = "ML")

Coefficients:
          ar1      ar2      ar3      ar4     sar1     sar2     sar3     sar4      xreg
      -0.8138  -0.4978  -0.4467  -0.3499  -1.4187  -1.2366  -0.8741  -0.4665  -19.4126
s.e.   0.0377   0.0449   0.0438   0.0380   0.0350   0.0571   0.0551   0.0340   52.4444

sigma^2 estimated as 195217:  log likelihood = -5327.59,  aic = 10673.17

BIC(model_ts1)

[1] 10720.78

Conclusion

From the above analysis we conclude that the model that better explains the above data is ARIMA.By changing the values (p,d,q) we can find the optimal value using the AIC and the BIC score which can help us in forecasting accurately.We learnt about dealing with trend and seasonality .The model can be improved further using residual components.

Team-9 members

Dona Sam -20BDA02

Pratibha KS - 20BDA15

Reba Susan- 20BDA37

Jayasree C-20BDA53

Abhijith - 20BDA60

Ananya Kumari- 20BDA68